library(party)
library(randomForest)
library(foreign)
## read "systemleveldata_r_t+1.csv"
nuke <-read.csv(file.choose())

## descriptive statistics
summary<-summary(nuke)
write.table(summary, file="summary.txt")

## conditional variable importance
data.controls <- cforest_unbiased(ntree=1000, mtry=2)

## To check the robusetness of results of conditional variable importance, the setting for "set.seed" is changes each time
set.seed(10507)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result1.txt")
abs(min(data.cforest.varimp))

set.seed(45)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result2.txt")
abs(min(data.cforest.varimp))

set.seed(9533)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result3.txt")
abs(min(data.cforest.varimp))

set.seed(6)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result4.txt")
abs(min(data.cforest.varimp))

set.seed(4276)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result5.txt")
abs(min(data.cforest.varimp))

## partial dependence plot
set.seed(10507)
aus.rf <- randomForest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, nuke, ntree=1000, mtry=2)
par(mfrow=c(3,2))
partialPlot(aus.rf, nuke, numnuke, xlab="# of nuclear states", main=NULL,  "versicolor")
partialPlot(aus.rf, nuke, nukeyear, xlab="Nuclear year counter", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, demonum, xlab="# of democratic states", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, gwp, xlab="Gross world product", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, unipolar, xlab="Unipolarity", main=NULL, "versicolor", axes=FALSE)
box()
axis(side = 1, at=(0:1))
axis(side = 2)
partialPlot(aus.rf, nuke, disstateratio, xlab="Lagged dispute-state ratio", main=NULL, "versicolor")

###############################################################

## fatal dispute-state ratio as the dependent variable
## conditional variable importance
data.controls <- cforest_unbiased(ntree=1000, mtry=2)

## To check the robusetness of results of conditional variable importance, the setting for "set.seed" is changes each time
set.seed(10507)
data.cforest <- cforest(fatalstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+fatalstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result1_fatal.txt")
abs(min(data.cforest.varimp))

set.seed(45)
data.cforest <- cforest(fatalstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+fatalstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result2_fatal.txt")
abs(min(data.cforest.varimp))

set.seed(9533)
data.cforest <- cforest(fatalstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+fatalstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result3_fatal.txt")
abs(min(data.cforest.varimp))

set.seed(6)
data.cforest <- cforest(fatalstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+fatalstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result4_fatal.txt")
abs(min(data.cforest.varimp))

set.seed(4276)
data.cforest <- cforest(fatalstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+fatalstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result5_fatal.txt")
abs(min(data.cforest.varimp))

## partial dependence plot
set.seed(10507)
aus.rf <- randomForest(fatalstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+fatalstateratio, nuke, ntree=1000, mtry=2)
par(mfrow=c(3,2))
partialPlot(aus.rf, nuke, numnuke, xlab="# of nuclear states", main=NULL,  "versicolor")
partialPlot(aus.rf, nuke, nukeyear, xlab="Nuclear year counter", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, demonum, xlab="# of democratic states", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, gwp, xlab="Gross world product", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, unipolar, xlab="Unipolarity", main=NULL, "versicolor", axes=FALSE)
box()
axis(side = 1, at=(0:1))
axis(side = 2)
partialPlot(aus.rf, nuke, fatalstateratio, xlab="Lagged fatal dispute-state ratio", main=NULL, "versicolor")

###############################################################

## war-state ratio as the dependent variable
## conditional variable importance
data.controls <- cforest_unbiased(ntree=1000, mtry=2)

## To check the robusetness of results of conditional variable importance, the setting for "set.seed" is changes each time
set.seed(10507)
data.cforest <- cforest(warstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+warstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result1_war.txt")
abs(min(data.cforest.varimp))

set.seed(45)
data.cforest <- cforest(warstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+warstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result2_war.txt")
abs(min(data.cforest.varimp))

set.seed(9533)
data.cforest <- cforest(warstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+warstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result3_war.txt")
abs(min(data.cforest.varimp))

set.seed(6)
data.cforest <- cforest(warstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+warstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result4_war.txt")
abs(min(data.cforest.varimp))

set.seed(4276)
data.cforest <- cforest(warstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+warstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result5_war.txt")
abs(min(data.cforest.varimp))

## partial dependence plot
set.seed(10507)
aus.rf <- randomForest(warstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+warstateratio, nuke, ntree=1000, mtry=2)
par(mfrow=c(3,2))
partialPlot(aus.rf, nuke, numnuke, xlab="# of nuclear states", main=NULL,  "versicolor")
partialPlot(aus.rf, nuke, nukeyear, xlab="Nuclear year counter", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, demonum, xlab="# of democratic states", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, gwp, xlab="Gross world product", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, unipolar, xlab="Unipolarity", main=NULL, "versicolor", axes=FALSE)
box()
axis(side = 1, at=(0:1))
axis(side = 2)
partialPlot(aus.rf, nuke, warstateratio, xlab="Lagged war-state ratio", main=NULL, "versicolor")

####################################################################

## robustness check - North Korea coded as a non-nuclear state
## conditional variable importance
data.controls <- cforest_unbiased(ntree=1000, mtry=2)

## To check the robusetness of results of conditional variable importance, the setting for "set.seed" is changes each time
set.seed(10507)
data.cforest <- cforest(disstateratiot1~numnuke2+nukeyear2+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result1_rob.txt")
abs(min(data.cforest.varimp))

set.seed(45)
data.cforest <- cforest(disstateratiot1~numnuke2+nukeyear2+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result2_rob.txt")
abs(min(data.cforest.varimp))

set.seed(9533)
data.cforest <- cforest(disstateratiot1~numnuke2+nukeyear2+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result3_rob.txt")
abs(min(data.cforest.varimp))

set.seed(6)
data.cforest <- cforest(disstateratiot1~numnuke2+nukeyear2+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result4_rob.txt")
abs(min(data.cforest.varimp))

set.seed(4276)
data.cforest <- cforest(disstateratiot1~numnuke2+nukeyear2+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result5_rob.txt")
abs(min(data.cforest.varimp))

## partial dependence plot
set.seed(10507)
aus.rf <- randomForest(disstateratiot1~numnuke2+nukeyear2+demonum+gwp+unipolar+disstateratio, nuke, ntree=1000, mtry=2)
par(mfrow=c(3,2))
partialPlot(aus.rf, nuke, numnuke2, xlab="# of nuclear states", main=NULL,  "versicolor")
partialPlot(aus.rf, nuke, nukeyear2, xlab="Nuclear year counter", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, demonum, xlab="# of democratic states", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, gwp, xlab="Gross world product", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, unipolar, xlab="Unipolarity", main=NULL, "versicolor", axes=FALSE)
box()
axis(side = 1, at=(0:1))
axis(side = 2)
partialPlot(aus.rf, nuke, disstateratio, xlab="Lagged dispute-state ratio", main=NULL, "versicolor")

###############################################################

## robustness checks by Christpher Way's (2012) coding of nuclear states
## conditional variable importance
data.controls <- cforest_unbiased(ntree=1000, mtry=2)

## To check the robusetness of results of conditional variable importance, the setting for "set.seed" is changes each time
set.seed(10507)
data.cforest <- cforest(disstateratiot1~numnuke_cway+nukeyear_cway+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result1_cway.txt")
abs(min(data.cforest.varimp))

set.seed(45)
data.cforest <- cforest(disstateratiot1~numnuke_cway+nukeyear_cway+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result2_cway.txt")
abs(min(data.cforest.varimp))

set.seed(9533)
data.cforest <- cforest(disstateratiot1~numnuke_cway+nukeyear_cway+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result3_cway.txt")
abs(min(data.cforest.varimp))

set.seed(6)
data.cforest <- cforest(disstateratiot1~numnuke_cway+nukeyear_cway+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result4_cway.txt")
abs(min(data.cforest.varimp))

set.seed(4276)
data.cforest <- cforest(disstateratiot1~numnuke_cway+nukeyear_cway+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result5_cway.txt")
abs(min(data.cforest.varimp))

## partial dependence plot
set.seed(10507)
aus.rf <- randomForest(disstateratiot1~numnuke_cway+nukeyear_cway+demonum+gwp+unipolar+disstateratio, nuke, ntree=1000, mtry=2)
par(mfrow=c(3,2))
partialPlot(aus.rf, nuke, numnuke_cway, xlab="# of nuclear states", main=NULL,  "versicolor")
partialPlot(aus.rf, nuke, nukeyear_cway, xlab="Nuclear year counter", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, demonum, xlab="# of democratic states", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, gwp, xlab="Gross world product", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, unipolar, xlab="Unipolarity", main=NULL, "versicolor", axes=FALSE)
box()
axis(side = 1, at=(0:1))
axis(side = 2)
partialPlot(aus.rf, nuke, disstateratio, xlab="Lagged dispute-state ratio", main=NULL, "versicolor")

####################################################################

## robustness check p(mtry)=3
## conditional variable importance
data.controls <- cforest_unbiased(ntree=1000, mtry=3)

## To check the robusetness of results of conditional variable importance, the setting for "set.seed" is changes each time
set.seed(10507)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result1_rob2.txt")
abs(min(data.cforest.varimp))

set.seed(45)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result2_rob2.txt")
abs(min(data.cforest.varimp))

set.seed(9533)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result3_rob2.txt")
abs(min(data.cforest.varimp))

set.seed(6)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result4_rob2.txt")
abs(min(data.cforest.varimp))

set.seed(4276)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result5_rob2.txt")
abs(min(data.cforest.varimp))

## partial dependence plot
set.seed(10507)
aus.rf <- randomForest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, nuke, ntree=1000, mtry=3)
par(mfrow=c(3,2))
partialPlot(aus.rf, nuke, numnuke, xlab="# of nuclear states", main=NULL,  "versicolor")
partialPlot(aus.rf, nuke, nukeyear, xlab="Nuclear year counter", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, demonum, xlab="# of democratic states", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, gwp, xlab="Gross world product", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, unipolar, xlab="Unipolarity", main=NULL, "versicolor", axes=FALSE)
box()
axis(side = 1, at=(0:1))
axis(side = 2)
partialPlot(aus.rf, nuke, disstateratio, xlab="Lagged dispute-state ratio", main=NULL, "versicolor")

###############################################################

## robustness check p(mtry)=4
## conditional variable importance
data.controls <- cforest_unbiased(ntree=1000, mtry=4)

## To check the robusetness of results of conditional variable importance, the setting for "set.seed" is changes each time
set.seed(10507)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result1_rob3.txt")
abs(min(data.cforest.varimp))

set.seed(45)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result2_rob3.txt")
abs(min(data.cforest.varimp))

set.seed(9533)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result3_rob3.txt")
abs(min(data.cforest.varimp))

set.seed(6)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result4_rob3.txt")
abs(min(data.cforest.varimp))

set.seed(4276)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result5_rob3.txt")
abs(min(data.cforest.varimp))

## partial dependence plot
set.seed(10507)
aus.rf <- randomForest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, nuke, ntree=1000, mtry=4)
par(mfrow=c(3,2))
partialPlot(aus.rf, nuke, numnuke, xlab="# of nuclear states", main=NULL,  "versicolor")
partialPlot(aus.rf, nuke, nukeyear, xlab="Nuclear year counter", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, demonum, xlab="# of democratic states", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, gwp, xlab="Gross world product", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, unipolar, xlab="Unipolarity", main=NULL, "versicolor", axes=FALSE)
box()
axis(side = 1, at=(0:1))
axis(side = 2)
partialPlot(aus.rf, nuke, disstateratio, xlab="Lagged dispute-state ratio", main=NULL, "versicolor")

###############################################################

## robustness check p(mtry)=5
## conditional variable importance
library(party)
data.controls <- cforest_unbiased(ntree=1000, mtry=5)

## To check the robusetness of results of conditional variable importance, the setting for "set.seed" is changes each time
set.seed(10507)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result1_rob4.txt")
abs(min(data.cforest.varimp))

set.seed(45)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result2_rob4.txt")
abs(min(data.cforest.varimp))

set.seed(9533)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result3_rob4.txt")
abs(min(data.cforest.varimp))

set.seed(6)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result4_rob4.txt")
abs(min(data.cforest.varimp))

set.seed(4276)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result5_rob4.txt")
abs(min(data.cforest.varimp))

## partial dependence plot
set.seed(10507)
aus.rf <- randomForest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, nuke, ntree=1000, mtry=5)
par(mfrow=c(3,2))
partialPlot(aus.rf, nuke, numnuke, xlab="# of nuclear states", main=NULL,  "versicolor")
partialPlot(aus.rf, nuke, nukeyear, xlab="Nuclear year counter", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, demonum, xlab="# of democratic states", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, gwp, xlab="Gross world product", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, unipolar, xlab="Unipolarity", main=NULL, "versicolor", axes=FALSE)
box()
axis(side = 1, at=(0:1))
axis(side = 2)
partialPlot(aus.rf, nuke, disstateratio, xlab="Lagged dispute-state ratio", main=NULL, "versicolor")

###############################################################

## robustness check p(mtry)=6, i.e. Bagging
## conditional variable importance
data.controls <- cforest_unbiased(ntree=1000, mtry=6)

## To check the robusetness of results of conditional variable importance, the setting for "set.seed" is changes each time
set.seed(10507)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result1_rob5.txt")
abs(min(data.cforest.varimp))

set.seed(45)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result2_rob5.txt")
abs(min(data.cforest.varimp))

set.seed(9533)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result3_rob5.txt")
abs(min(data.cforest.varimp))

set.seed(6)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result4_rob5.txt")
abs(min(data.cforest.varimp))

set.seed(4276)
data.cforest <- cforest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, data = nuke, controls=data.controls)
data.cforest.varimp <- varimp(data.cforest, conditional = TRUE)
data.cforest.varimp
write.table(data.cforest.varimp, file="result5_rob5.txt")
abs(min(data.cforest.varimp))

## partial dependence plot
set.seed(10507)
aus.rf <- randomForest(disstateratiot1~numnuke+nukeyear+demonum+gwp+unipolar+disstateratio, nuke, ntree=1000, mtry=6)
par(mfrow=c(3,2))
partialPlot(aus.rf, nuke, numnuke, xlab="# of nuclear states", main=NULL,  "versicolor")
partialPlot(aus.rf, nuke, nukeyear, xlab="Nuclear year counter", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, demonum, xlab="# of democratic states", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, gwp, xlab="Gross world product", main=NULL, "versicolor")
partialPlot(aus.rf, nuke, unipolar, xlab="Unipolarity", main=NULL, "versicolor", axes=FALSE)
box()
axis(side = 1, at=(0:1))
axis(side = 2)
partialPlot(aus.rf, nuke, disstateratio, xlab="Lagged dispute-state ratio", main=NULL, "versicolor")

